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Objective 

Currently no software simulations are available from industry or 
academia which can allow an accurate computation of the radar signature 
from full scale non-metallic vehicles primarily because exl ^ 1I |S 
formulations have excessive memory and computational demands 
addition to being unsuitable for modem machine architectures Such 
software are required for test and evaluation purposes of future aircraft 
systems/platforms, target simulation platforms, on-board communication 
systems and for radar signature evaluation. 

Traditional methodologies for radar signature or performance 
analysis of on-board communication systems have not shown any promise 
^simulating large scale composite vehicles. Furthermore, approximate 
simidadon technfques are of no utility in modeling complex compose 
structures. The intent of this research is thus aimed at developing 
methodologies which will allow accurate and efficient simulations of ful 
scale vehicles. We propose to accomplish this with *e developmem of 
algorithms which have very low memory requirements and ^make use of 
recent advances in high performance and parallel f . 

Darticular technique to be employed is a hybrid version of the fin 
element method which has been successfully implemented and tested at the 
University of Michigan for two-dimensional applications. Als0 ’ s0 ™ ‘ 
dimensional results relating to specific scattermg and antenna 
configurations have verified our expectations of the 
The basis of this approach is to reformulate an open-boundary p 
into a closed-boundary one using a mesh termination method which b g 
the outer boundary very close to the body. Segmentation and c °nnecuv ty 
schemes can then be employed to subdivide the system/volume 
(of any number) sections which can be independently treated by separate 
single o? parallel processors. Except for the small Toepl.tz connect, v.ty 
matrices, all others are sparse and banded, making their solution 
convenient, and since there is no need for concurrent storage of any of 
these submatrices, there is no limitation on the size of the computationa 
domain, a problem which has become the bottleneck in computational 
electromagnetics. It is thus our claim that this solution approach, coupled 
withT aStages of the finite element method, will make a dramatic 
impact in the field of numerical electromagnetics, and promises to provide 
the necessary future tools for numerous applications, including target 
stauladon radar signature analysis and imaging, antenna modeling and 
VLSI circuit simulations. 
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Progress 


As stated above, the main goal of this project is to develop methodologies 
for scattering by airborne composite vehicles. Although our primary focus 
continues to be the development of a general purpose code for analyzing 
the entire structure as a single unit, a number of other tasks are also 
pursued in parallel with this effort. These tasks, are important in testing 
the overall approach and in developing suitable models for materials 
coatings, junctions and, more generally, in assessing the effectiveness of the 
various parts comprising the final code. 

Below we briefly discuss our progress on the five different tasks which 
were pursued during this period. Our progress on each of these tasks is 
described in the four detailed reports (025921 -24-T, 025921 -25-T, 
025921 -26-T, 025921 -27 -T) and the memoranda included with this 
document. The first task described below is, of course, the core of this 
project and deals with the development of the overall code. Undoubtedly, 
it is the outcome of the research which was funded by NASA-Ames and the 
Navy over the past three years. 


Task: Three-dimensional hybrid finite element formulation for 
scattering 


During this period we reached an important milestone in the development 
of the proposed electromagnetic tools. Specifically a full three dimensional 
code has been written and validated for computing the scattering by open 
bodies having arbitrary shape and material parameters. The completion of 
this code required the successful development of several new theoretical 
and mathematical procedures. Most importantly, a new mesh termination 
scheme was established to ensure the outgoing nature of the waves. With 
this scheme, the termination boundary was placed very close to the target 
thus greatly minimizing the unknown count. Moreover, the interior region 
was discretized with a new set of edge-based finite elements which do not 
suffer from the instabilities of the traditional node-based elements. By 
virtue of the employed near neighbor finite element method, the resulting 
system is sparse and when solved via an iterative method it has minimal 
memory requirements. We have so far solved systems having an unknown 
count of over 200,000 and our goal in the next few months is to increase 
the system size to more than 500,000. A major effort has also be devoted 
in establishing pre and post processing software and interfacing these with 
the new analysis code. As of the moment, we have interfaced our codes we 
a major commercial geometry/mesh generation package and several color 
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graphics tools and GUI platforms for entering and displaying the input and 
output data. More information on this task are provided in the University 
of Michigan technical report 0259121-24-T entitled "Application of Edge- 
Based Finite Elements and Vector ABCs in 3D Scattering. Representative 
results extracted from this report are shown in figures 1 and 2. 1 hese 

correspond to a metallic cube and circular cylinder, respectively. 


Task: Improvements of the hybrid finite element method for a 
body of revolution 


Last year a new finite element-boundary integral formulation for a body of 
revolution (BOR) was developed and coded as described in the University 
of Michigan technical report 025921-1 8-T. The developed formulation 
was superior from others in being capable of analyzing bodies with 
material coatings and inhomogeneous sections. In addition, by virtue ot the 
formulation, only O(n) storage was required, important for large scale 
computations. However, it soon became apparent that the developed code 
was unsuitable for large bodies because of resonance conditions associated 
with the discrete system. The resonance characteristics are due to the 
boundary integral and have been a perennial difficulty to all formulations 
which make use of this integral. During this period we investigated 
several new method for eliminating these resonances and improving the 
efficiency of the code. 


First, a new absorber termination (ATB) scheme was developed to replace 
the boundary integral (see figure 3). This approach leads to fully sparse 
matrices and is inherently rid of the resonances. However, the ATB had to 
be placed more than 1/2 wavelengths away from the structure, thus 
compromising the efficiency of the solution. At the moment we are 
investigating improvements to the ATB to be reported during the nex 

period. 

A second method for eliminating the resonances involved a recasting of the 
boundary integral with the scattered (rather than the total) field as the 
unknown. This approach is quite successful as illustrated in figure 4 tor 
eliminating the resonances in two dimensions and is currently being 
implemented for three dimensional computations. 

A third approach based on the implementation of a "complexification" 
scheme was also found equally successful in eliminating the resonances. 
This approach is attractive because of its simplicity. More details on its 
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implementation will be reported during the next period along with results 
and a description of the final code. 


Additional details on our progress pertaining to this task are given in the 
University of Michigan technical report 02592 1-26-T. 


Task: Resistive card modeling in the context of finite element 

solutions 

Resistive cards are nowadays used on many airborne vehicles for radar 
cross section control. It is thus important to provide appropriate models 
for resistive cards in the context of the finite element method which is the 
basis of our methodology. The difficulty with traditional modeling of 
resistive cards is the implementation of the jump condition on the field 
across the sheet. This requires the specification of double nodes at the 
location of each resistive sheet, resulting in cumbersome and somewhat 
inefficient formulations. To avoid the double node specification, we 
modified the functional which is normally minimized in the implementation 
of the finite element method. The new functional does not require the 
specification of double nodes at the location of the resistive sheet, making it 
attractive to use in conjunction with commercial mesh generation packages. 
The new functional was tested and implemented in the new code described 
in the first task above. Details of the formulation are discussed in the 
University of Michigan report 025921 -27-T. 


Task: New finite element formulations and mesh termination 
schemes 

In parallel with the development of the above more traditional hybrid finite 
element formulation and associated code we are also exploring new 
methods of implementation. The goal with these new implementations is 
efficiency, accuracy and simplicity. For example, in order to increase the 
efficiency of the solution we have considered new expansion functions and 
shapes of elements with the intent of using different expansion functions 
within a single mesh. The idea of this implementation is to provide the 
formulation with as much apriori information as possible, thus reducing 
the computational intensity and unknown count without comprising 
accuracy and geometrical adaptability. For this purpose, a two dimensional 
code was written which may break up the computational domain in 
different regions depending on the expected field complexity in those 
regions. Those regions over flat surfaces which are dominated by 
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geometrical optics fields can be modeled with large size elements and field 
shaDe functions which include in their representat.on the known 
geometrical optics o PTD fields. Small correction «ms can *en be 
introduced to provide the complete expression for the field in those 

regions. 

Furthermore to remove the requirement of a boundary integral for 
terminating the mesh we have separately explored the possibility of 
a“ive mesh terminations schemes. The employed ide* . relies on the 
use of artificial absorbers placed in the near proximity of the body^ The 
advantage of these absorbers is that they retain the sparsity of the matrix 
and thif is crucial in reducing the computational intensity and memory 
requirement of the method. We have already developed and tested such 
absorbers (ATBs) and we continue to improve their absorptive properties 
to bring *em as close to the body as possible without comprising accuracy. 

Some of the details pertaining to this task are given in the included 
memorandum by D.C. Ross and J.L. Volakis dated January 5, 1992. 

Task: Applications of wavelets to electromagnetics 

Wavelets are a family of basis functions which are generated by translation 
and dilation of a single function. In the past they have been used for 
compressing data or functions defined over a large range. Hus is a° n e by 
weTghSntegratmg or otherwise transfoiming the data wit . each o the 
basis in the wavelet family, thus generating an new array of data " h 
referred as the wavelet transform of the ongmal data or funcuon. The 
attractive property of the wavelet transform is its bandlimited form. This 
"ause through inverse transformation each of the wavelet basis or 
functions adds derails to the original function that are ^“"" a ^ a °^ 
in accordance with their level of importance. More on the fundamenta 
of this property is described in the attached memorandum y • P * 

In this P memorandum an application of wavelets to f 

described. In particular the wavelet transform is instead applied to a large 
complex matrix resulting from the discretization of a typical boundary 
integral. By application of the wavelet transform, it is shown that a full 
matrix can be compressed to a banded matrix, important for reducing t 
solution's memory requirements. The included memorandum is a 
preliminary study and additional applications are planned m the next te 
months The goal is to take advantage of the wavelet transform in treating 
the boundary integral matrix arising in some hybrid finite element 

formulations. 


5 


Task: Alternative integral equations for scattering by three 
dimensional coated structures 


Thin material coatings are nowadays routinely used in the construction all 
airborne vehicles. Nevertheless, their electromagnetic characterization in 
the context of integral equations still relies on the use of the simple but less 
accurate impedance boundary conditions. On the other hand, although 
finite element formulations can readily deal with the presence of material 
coatings, their implementation is inefficient when the coating becomes thin. 

This is because very small elements must be used in the vicinity of the thin 
coating. 


To avoid the inaccuracy of the impedance boundary condition and the 
inefficiency of the finite element method for modeling thin coatings, a new 
integral equation was derived and discretized for implementation. The 
derivation of the new integral equation is described in the University of 
Michigan technical report 02592 1-25-T entitled "Alternative field 
representations and integral equations for modeling inhomogeneous 
dielectrics. In comparison with traditional integral equations which 
require as many as eight unknowns in their discretization, the new integral 
equations requires only three unknowns per cell without any compromise 
m accuracy. The approximate integral equation based on the impedance 
boundary condition requires four unknowns per cell. In addition the new 
integral equation is also associated with a lower kernel singularity which 
greatly facilitates its discretization. A code based on this formulation is 
currently being written and results based on this code will be reported 
during the next few months. 


Task. New impedance wedge diffraction coefficients 

Sometime ago we reported in the University of Michigan report 025921-4- 
T a new approximate impedance wedge diffraction coefficient for 
characterizing the diffraction by an impedance wedge. We resorted to an 
approximate coefficient only because the impedance wedge is not amenable 
to an exact solution when excited by a wave incident in a plane other than 
that normal to the wedge. This diffraction coefficient was later 
implemented and successfully tested in a traditional PTD code under a 
contract supported by the U.S. Army (see University of Michigan report 
028371-1-F). 
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After several applications, the ^rosl'polari^onfields. 

found to be at times inaccurate for ^T?V E ^ ibi i it y of a more rigorous 
To remedy this deficiency, we From the start, this 

development for the proposed dlffrac «°" ^ arrived at a procedure which 
effort proved challenging u coefficient The new coefficient is 

has yielded a more accurate faction upcoming technical report. 

S'SSdSa 1 - — d '- ”■ 

currents for possible implementation in a PTD code. 


Transitions 

During the next period we shall emphasize die following 

1 . -nte application of the 

S"sSTprov" the method's suitability for large scale 
implementations. 

2. Develop a segmentation and/or 

ap^oach'shouk^focilitate the discretization of large size structures and 
computations associated with full scale vehicles. 

*■ sr 

thus* increase the computational efficiency of the method. 

5 - 

difficulties associated with the internal body resonances. 

6. Complete and validate the code based the new integral equation for 
coated bodies. 

7 Complete the development and test the accuracy of the new skew 
incidence impedance wedge diffraction coefficient. 


7 


List of reports submitted this period 

A. Chatterjee, J.M. Jin, and J.L. Volakis, "Application of edge-based finite 

e ements and vector ABCs in 3D scattering," Univ. of Michigan Radiation 
Laboratory Report #02592 1-24-T, Jan. 1992. 8 Kaaiation 

mnH^i° IaklS ’u Alternative field representations and integral equations for 
modeling inhomogeneous dielectrics," Univ. of Michigan Radiation 
Laboratory Report #025921 -25-T, Feb. 1992. 8 Radiatl0n 

fo^'scaneirhv h J 0laki f S ’ " Pr °, gKSS °" hybrid finite element methods 
scattering by bodies of revolution," Univ. of Michigan Radial,, 

Laboratory Report #02592 1-26-T, Feb. 1992. h 8 Radiation 

X° lakis ’ C L - Yu and A C Wo °. "Modeling of resistive 

U " iv of 1 Michigan Radiati °" Lab -^ 

Semiannu a| Progress Report for NASA Grant NAG-2-541 "Develonment 
of 3D Electromagnetic Modeling Tools for Airborne Vehicles ™Unh? of 
Michigan Radiation Laboratory Report #025921-28-T, Feb. 1992 


8 



Bistatic pattern of a .755 X metallic cube 


0 =180" 






(£ =90* cut 
• Measured 
cut 

\* Measured 


60 90 120 

Observation Angie 0 o , deg. 


to _ iso 0 ! of a metallic cube having an edge length 
Fi6U ' e 11 2»d order vector absorbing boundary condition. 




Backscatter echo-area/A*. dB 



Backscatter pattern of pec cylinder (r=3A,h=6A) 



Figure 2: 


Backscatter pattern of the cylinder in Figure 7 enclos. 

mal outer boundary and using the second order absoi 
condition. 


a confor- 
boundary 





Echo Width [dB| Echo Width [dB] 



(a) 



(b) 


Figure 3: Backscattering echowidth of a circular cylinder: (a) original formula 
tion; (b) scattered field formulation with complexification. 
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A Finite Element Formulation for Scattering from 
Electrically Large 2-Dimensional Structures. 


Daniel C. Ross and John L. Volakis 


January 5, 1992 


Abstract — A Finite Element Formulation is given using the scattered field ap- 
proach with a fictitious material absorber to truncate the mesh. The formulation in- 
cludes the use of arbitrary approximation functions so that more accurate results can 
be achieved without any modification to the software. Additionally, non-polynomial 
approximation functions can be used, including complex approximation functions. 
The banded system that results is solved with an efficient sparse/banded iterative 
scheme and as a consequence, large structures can be analyzed. Results are given 
for simple cases to verify the formulation and also for large, complex geometries. 


I. Introduction 


The popularity of the Finite Element Method (FEM) in Electromagnetic Scattering is largely 
due to the fact that the method results in a banded system with an O(N) memory demand. This 
is in sharp contrast to the Moment Method which has an 0(N 2 ) memory demand. For open 
region scattering problems however, the Finite Element region must be artificially truncated in 
order to bound the problem. This has been accomplished with approximate methods such as 
the Uni-Moment Method [1] and the Absorbing Boundary Condition [6] and more successfully 
with the rigorous Boundary Integral Technique [3]. Unfortunately, the solution of the Boundary 
Integral in general results in a dense system thus destroying the O(N) character of the FEM. 
Also, the boundary integral technique breaks down at resonant frequencies and these resonant 
frequencies become more and more closely spaced as the size of the boundary is increased. 

More recently, the technique of using a fictitious material absorber to confine the scattered 
field has been suggested in [4]. Since fictitious materials having a negative real part of n and 
c can be used, an absorber with a very low reflection coefficient even at near grazing incidence 
can be realized. Also, this technique requires the unknown quantity to be the scattered field, 
not the total field as is usually the case. The benefits of this are two-fold. First, the scattered 
field does not vary as rapidly as the total field and second, the scattered field approach requires 
that the incident field be evaluated in the interior of the mesh, thus reducing the phase error 
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Figure 1: Region of solution 

normally associated with the total field approach. That is, the total field approach requires that 
the incident field be evaluated only at the fictitious boundary and then propagate through the 
mesh to the scatterer thus introducing phase error, especially notable when the mesh is large. It 
has been confirmed in [5] that the scattered field approach results in a more accurate solution. 

It is known that there is also an overall inherent error in the FEM that does not exist in the 
Moment Method [7], This error is a result of the bandedness of the system and is introduced 
as the solution propagates through the mesh. In fact, experience has shown that as the mesh 
area increases, the element size must decrease. The only solution to this problem is to use higher 
order approximation functions, which usually involves a complete re-formulation and re-design 
of the software. 

In this paper a number of advances in FEM calculations are combined to give a solution 
to the 2-Dimensional Helmholtz equation for electrically large scatterers. The scattered field 
approach is used and a fictitious material absorber is used to truncate the mesh. Unique to this 
formulation is the use of arbitrary approximation functions that make it possible to increase the 
accuracy of the solution without modifying the software. 


II. Formulation 


Consider the region shown in Figure 1. 

Given that there is some incident field, we wish to solve the wave equation in the region $7 
and find the fields scattered by the dielectric/magnetic/conducting target. 
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( 1 ) 


The source free wave equation in two dimensions can be written 

V ■ (uV<}>)+ klv<t> = 0 


where 


u = 1//J r 

v = l/e r 

for E polarization ( z directed electric field) 


and 


<t> =■ H z 

U = l/f r 

V = 1/flr 

for II polarization ( z directed magnetic field) 


( 2 ) 

(3) 

(4) 


(5) 

( 6 ) 
(7) 


Now let T be some compactly supported testing function used to enforce (1) over the entire 
region. 



r[V • (i/V0) 


4- ds — 0 


( 8 ) 


Using the identity 

T[V ■ (t/V^>)] = ■ VT) + V • ( TuV4> ) (9) 


gives 



[u(vr- v<f>) 


k1vT4>] ds 



V • ( TuV<t>)ds 


= 0 


( 10 ) 


and applying the divergence theorem to the second term in the above brings us to the weak form 
of eqn . ( 1 ) . 



[u(VT • V<j>) 



uT^dt- (f uT^-di = 0 

dn J Tc on 


(H) 


To use eqn. (11) in an FEM implementation would be considered a total field approach since 
the unknowns would be the total field. To get the analog of (11) with the scattered field as the 
unknown we let 


<j> = 4> s + <t>' 


( 12 ) 


3 


and (11) becomes 


j JnVT V <!>’)- klvT6’]ds - j f T(u^ + u^-)de (13) 

-jf T(u ^ + i ^£ )de = - J J n MVT -V#) - klvTtfds 

where u and v are the material constants of the scatterer in a free space environment while 
u and v include the fictitious absorber. In this way the incident field does not interact with the 
fictitious absorber. 


In either case, eqns. (11) and (13) require that the normal derivative of the unknown (total 
field and scattered field respectively) be specified on the enclosing boundaries r a and T c . These 
normal derivatives are arbitrary, and once given, a solution in the interior is uniquely defined. 

For E polarization we will put a perfect magnetic conductor at r o and for H polarization we 
will put a perfect electric conductor at T a thus the second term on the LHS of (13) is 

-l* T t dl <“> 

Note that this conductor is seen by the scattered field only, not the total field. 


The boundary condition at r c is 


d<t>’ dp 

^r + ^r = ° 

for II polarization, and thus the third term on the LHS of (13) is zero. 


(15) 


For E polarization the boundary condition at T c is 

P + P = o (i6) 

and since the solution is specified on this contour, the normal derivatives will uncouple from the 
rest of the equation once we have discritized (13) and created a linear system, thus the third 
term on the LIIS of (13) drops out. 


We have thateqn. (13) becomes 

[ [ [l/Pr(VT -VP)- kl( r Tp) ds - / l/^rffr 
J •'O Jp m on 

- yV[l/£r(Vr . Vp) - kprTp] ds 

for E polarization and 

/ / [l/c r (vr • V^) - klurTp) ds - f 1 !i r T d 4-dt = 

J •'H J p a on 

~ ! ln [l,lr[VT ' V0,) “ ds 


(17) 


(18) 
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for H polarization. 

The region ft now is divided into arbitrary, polygonal finite elements. Each element is defined 
by its corner nodes and may contain other nodes located anywhere in the element. The solution 
in the e*th element is approximated by 


N* 


4> c (x,y) = ^TS'(x,y)«£ 


(19) 


n = 1 


where N' is the total number of nodes in the e’th element, 4>' n is the value of * at the n th node 
and 55 (x, y) is the n’th shape function for the e’th element understood to be zero outs.de of lie 
e’th element. These shape functions will be defined in the next section. To solve (17) and (18) 
using Galerkin’s method, we let the testing function T be the same as the basis function in (19). 


Substituting (19) into (17) and (18) gives the Element Equations 

N‘ 


IS, 


( 20 ) 


( 21 ) 


[i/^vs, e ■ v Y S'4>] - kit' St Y s]4>j] ds = 

e j=l j=l 

+ ^ 1/tiS'^dC- Jh [l/X(VS! -VW-kteS'Wds i= 1,2,3...* 

for E polarization and 

[ \ [l/e*VS* • V Y S i*i ~ k ° Y % ds = 

J j — 1 j — 1 

+ j> l/c'S'^£de- J [1/£*(V5/ -VW-klftS'^ds i= 1,2,3...* 

for H polarization. 

Again, ft and l\ are the material constants of the scatterer in a free space environment while 
and c e r include the fictitious absorber. 

Equations (20) and (21) define an N' by N e linear system for a single element. Assembling 
all of the element systems together and converting to global node numbers gives the global finite 
element system 


(*]{**} = {*■} 


( 22 ) 


For E polarization only, the boundary condition (16) must be enforced by setting 4>‘ - -<f>' 
for nodes on the conductor. 


III. Arbitrary Approximation Functions 
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In eqn. (19) the expansion of the solution was given as a sum of A 7 shape functions weighted 
by the nodal values of the solution. To determine these shape functions let 


N e 

<t> e (x,y) = ^a n / n (x,y) ( 23 ) 

r» = l 

where the / n (x,t/)’s are linearly independent functions known as approximation functions. 


Forcing <f) e (x,y) to be <t>*(Xi,yi) at each nodal point gives the system 

W5JW} = m ( 24 ) 

where [F t *] is an N by N matrix whose i,j y th element is given by /;(x,-,t/,) and x, , y t are 
the coordinates of the Fth node. For ease of programming, these approximation functions are 
often chosen to be linear, ie. f x = l,/ 2 = x,/ 3 = y. However, since it is known that higher 
order approximation functions result in better conditioned systems and are more accurate 1 , it is 
desirable to use approximation functions of arbitrary order or even of a non-standard type such 
as complex transcendental functions which are better suited for Electromagnetics. 

The shape functions are then given by 


N* 

Sn(X’y) = '52fdx,y)[F'„]- 1 (25) 

i=l 

Notice that the first derivatives of the shape functions are also required and this must be taken 
into account when choosing approximation functions since from (23), derivatives of the shape 
function transfer directly to the approximation functions. 


Note also that the inverse of an N e by N e matrix is required to calculate the shape function 
for each element. Since N e will not be very large, this does not introduce any great expense. 

When the arbitrary shape functions are used to generate the element matrices an efficient 
numerical integration scheme must be used. Since the elements are arbitrary polygons, an in- 
tegration algorithm based on the polygon fill algorithm from computer graphics was developed. 
This integration algorithm takes as input: the number of points defining the polygon, the coordi- 
nates of these points, the number of sub-regions to split the area into, the order of the Gaussian 
quadrature to use, and the integrand function. With some experimentation, it was found that 
increasing the accuracy of the integrations had little effect on the far field scattering patterns 
once the integrations were approximately five percent accurate (for linear shape functions.) 


IV. Results 


The bistatic RCS pattern for a perfectly conducting circular cylinder of radius A/2 is shown in 

r won' Lmear ShapC functions were used and “ can be seen from figure 3 a discretization 
of A/20 was required. The fictitious absorber was placed A/2 away from the cylinder to avoid 
interactions with the scatterer. 
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Figure 4: Bistatic RCS of Airfoil 


A more useful result is shown in figure 4 for a perfectly conducting airfoil at low frequency. 
The wing is approximately l|A long and as can be seen from figure 5 a discretization of A/20 
w*as again required but the fictitious absorber needed to be placed a full wavelength away from 
the target in order to get accurate results for H polarization. The reason for the discrepancy in 
the E polarization results is that the artificial absorber actually needs to be placed farther away 
from the scatterer. The following discussion addresses this issue. 

For quadratic shape functions, which are obtained by setting f\ = l,/ 2 = x , / 3 = y, f 4 = 
x y*/s = z 2 , and / 6 = y 2 , the mesh is not as dense as can be seen from figure 6. However since 
there are twice as many nodes per element, the net effect is that there are about half as many 
elements but the same number of nodes. However the increase in accuracy makes it possible to 
analyze larger and more detailed scatterers. 

For cubic shape functions, which are obtained by setting, for example, f x = l,/ 2 = x ,/ 3 = 

V,U = xy>h = = y 2 ,/? = = y^/s = the mesh is even less dense as can be 

seen from figure 7. Again, there are approximately the same number of nodes but the number 
of elements is cut in half. Actually, since the accuracy is increased exponentially, the mesh size 
could be increased further thus actually reducing the number of unknowns. This would seem 
to imply that higher order approximation functions keep getting more efficient. Whether or not 
this is the case remains to be seen. 


Another choice of approximation function could be the Ilankel function of order 0. This would 
make it possible to discretize using larger elements far away from the scatterer where the wave 
becomes somewhat cylinderical. There would still be other approximation functions to allow 

1 Polynomial approximation functions of Order p result in a solution which is h p accurate where h is the mesh 
size. 
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Figure 5: Mesh for Airfoil : Linear 



Figure 6: Mesh for Airfoil : Parabolic 
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Figure 7: Mesh for Airfoil : Cubic 
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Abstract 

Wavelets are an exciting new topic in applied mathematics and 
signal processing. This paper will provide a brief review of wavelets 
which are aiso known as families of functions with an emphasis on 
interpretation rather than rigor. We will devise an indirect use of 
wavelets for the solution of integral equations based upon techniques 
adapted from image processing. Examples for resistive strips will be 
given illustrating the effect of these techniques as well as there promise 
in reducing dramatically the storage requirement in order to solve an 
integral equation for large bodies. We also will present a direct im- 
plementation of wavelets to solve an integral equation. Both methods 
suggest future research topics and may hold promise for a variety of 
uses in computational electromagnetics. 


1 Introduction 

Wavelets have generated significant excitement among applied mathemati- 
cians and engineers recently due to their unique properties and diversity. 
Wavelets are actually a family of basis functions which are generated by 
translation and dilation of a single function. The uses are diverse since many 
different wavelets may be generated which possess some properties which 
are particularly well suited for a given purpose. Three prominent uses of 
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wavelets are image compression which is discussed by Mallat ( [1 ] , [2] ) and 
Daubechies[3], solution of partial differential equations by the staff at Aware 
Inc. ([4], [5]), and the solution of Fredholm Integral Equations by several 
members of the Yale research staff([6],[7],[8]). 

This paper will give a brief and basic review of wavelets with an empha- 
sis on the character of these functions rather than a rigorous mathematical 
treatment. Our purpose is to apply wavelet expansions or wavelet trans- 
forms to the solution of electromagnetic problems. We will show that a novel 
technique which is adopted from image processing may be applied to tra- 
ditional Method of Moments(MoM) computer codes to achieve significant 
storage savings. Since modern computers are becoming increasingly fast, the 
primary limitation on MoM implementations is the 0{N 2 ) storage require- 
ment. The use of wavelets is intended to alleviate this requirement. This 
“image processing” technique may be used with existing MoM implementa- 
tions to reduce storage. Finally we will look at direct applications of wavelets 
in electromagnetics with a modal solution of the Fredholm equation and dis- 
cuss future research topics. Although we have only used E-polarized resistive 
strip examples for this paper these techniques are applicable to other elec- 
tromagnetic problems of interest. In particular, although we only use MoM 
examples herein, we are interested in applying wavelet analysis in the solution 
of the boundary integral associated with the FEM-BE method. 


2 Wavelets and Dilation Equations 

Wavelets and their underlying form, a dilation equation, is an unfamiliar 
concept to most engineers. We will attempt to give a feel for what a wavelet 
is and what properties it possesses which are attractive for electromagnetics 
without going into detail. Strang[9] gave an excellent high-level review of 
the properties of wavelets while Daubechies[lO] provides a rigorous review 
with significant new material. Mallat[2] was able to relate wavelet decompo- 
sition and reconstruction in terms of multiresolution analysis which allowed 
significant progress to be made in exploiting wavelets for the aforementioned 
applications. 

In order to understand wavelets, we must first look at a dilation equation 

H x ) = c k<i>{2x - k) ( 1 ) 
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which is recognized as a two-scale difference equation. Strang[9] reviews 
the criterion which must be met in order for ( 1 ) to be unique. The “scal- 
ing function", <f>(x), is determined recursively and thus many properties of 
it must be observed rather than defined. Daubechies[lO] was able to con- 
struct orthonormal wavelets with compact support in both domains which 
was thought to be impossible by using the recursion 

< M* r ) = XI (~ X _ ^ ') 

k 

with the box function as the fundamental function (©oU)) and certain simple 
conditions which determine c*. The family of wavelets derived by Daubechies 
will be used for the examples of this paper although this is certainly not the 
only or necessarily optimal choice of wavelets. 

At this point, we introduce the continuous wavelet family which comes 

from the scaling function, 0(x), as 

i'(x) = ^4>{s{x - u)) (3) 

where the scale parameter(s) and the translation parameter(u) range over the 
positive real axis and the entire real axis, respectively (The notation adopted 
for the remainder of this paper is from Mallat [1]). A fundamental property 
of wavelets is their approximation capability. Daubechies family is denoted 
by the number of recursions used to determine the analyzing wavelet. For 
example we will denote the Daubechies wavelet which has four coefficients 
as DAUB4 which can recover up to a linear polynomial since the two lowest- 
order moment vanishes as discussed by Strang[9]. As terms are increased, 
more moments vanish and thus the polynomial approximation increases in 
order (e.g. DAUB6 will recover a quadratic). 

As was previously mentioned, the Daubechies family is not only orthonor- 
mal but they enjoy the useful property of compact support. For example, 
figure 1 illustrates DAUB20 wavelet compared with the comparable cosine 
if we wish to represent an impulse “frequency" of / = 10 24 A where A is the 
spatial sample size. As seen, the cosine(Fourier series basis) is compact in 
the “frequency” domain yet requires infinite support in the spatial domain. 
However, the Daubechies wavelet is compact in both domains. A further 
example is given in figure 2 where the excitation impulse is at / — fo 24 A 
and we observe that as the “frequency" increases, the analyzing wavelet has 
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diminishing spatial support. This is the most important property of wavelets 
since finer spatial resolution is used for “high frequency” components which 
results in a more accurate evaluation than is possible for the comparable 
Fourier coefficient which is computed with constant sample size regardless of 
the “frequency”. 

Several groups ([4], [5], [6], [7], [8]) have recognized that the property of 
compact support in both domains will result in a sparse matrix when an 
integral operator is converted to a matrix operator. We will discuss in the 
next section the fact that if one thinks in terms of multiresolution analysis, 
the preceding observation is obvious. 


3 Matrix Compression Techniques 

One of the most promising applications of wavelets is image compression 
for the transmission and storage of detailed images. For example, The FBI 
is looking at using wavelet compression techniques in order to store their 
extensive finger print library without loss of fidelity. We will give a brief 
review of the wavelet transform, multiresolution analysis , and what this 
topic has to do with electromagnetics. 

The continuous wavelet transform is given here without proof[l] 

F(s,u) = J^f(x)xf) t (x-u)dx (4) 

where the wavelet family is given 

&(*) = y/s^j{sx) ( 5 ) 

and 4>(x) is the mother wavelet which is denoted as <f>{x) in (8). In order to 
reconstruct a signal from its wavelet transform, we define the inverse wavelet 
transform as 


f( x ) - — u)dsdu ( 6 ) 

A more practical transform for numerical applications is the Discrete Wavelet 
Transform (DWT) which acts on a sequence of samples. In order to cover the 
phase-space, we need to uniformly sample the translation parameter while 
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exponentially sampling the scale parameter! 1]. The resulting DWT is given 

Wi>[/(j,n)] = / f(x)v 21 {u - rm~ J )dx (7) 

where a is the elementary dilation step. The signal may he reconstructed 
from its DWT by the series 

/( i )=mz < f{u),4'v{ u - > v 2>i u - nQ ~ 3 ) ( § ) 

j n 

with < > denoting the standard inner product. In practice, we choose 

a = 2 for historical reasons although this is not necessary. 

The Fourier transform is recognized as a transformation from the spatial 
domain to the wavenumber domain. However, the wavelet transform is a 
transformation from the spatial domain to the scale domain. This is the heart 
of the multiresolution interpretation of (7) as well as its fast implementation, 
the Fast Wavelet Transform(FWT). Each successive coefficient of the DWT 
describes the difference or details between the current level and a the previous 
decimated level. Figure 3 illustrates this concept. The first level, figure 
3a, describes the smooth component of the signal. Each successive level 
adds detail to the reconstruction. Thus a smooth signal will have a large 
large-scale component and very small detail components which results in 
signal compaction. The efficiency of the FW r T algorithm follows from the 
pyramid structure shown in figure 3. The highest detailed level (smallest- 
scale) requires that all samples be used in its computation. However each 
successive level requires half (a = 2) the number of samples due to dilation 
which results in the FWT having the same order of operations as the FFT. 

Mallat[2] was the first to recognize that this property naturally satis- 
fies the procedure used in pyramid-encoding schemes. If one applies a two- 
dimensional FWT which is computed in a manner similar to the most com- 
mon implementation of the two-dimensional FFT to an image, the image 
can be compressed into few large components and many components which 
are less than some tolerance whose omission does not adversely affect the 
fidelity of the image. Significant compression can occur depending on which 
analyzing wavelet is used. 

We realized that the typical impedance matrix generated by a MoM code 
may be thought of as a complex- valued image! Thus it is reasonable to as- 
sume that we may apply the wavelet transform to a dense impedance matrix 
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to generate a sparse matrix. It should be mentioned that this technique 
is similar in concept to the Impedance Matrix Localization(IML) method 
proposed by Canning[ll] who uses physical arguments to generate a trans- 
formation matrix which yields a sparse matrix by using sub-domain basis 
functions. It is not yet apparent which method yield superior results. 

Suppose one wants to solve for the E-polarized scattering from a resistive 
strip. The appropriate integral equation is a Fredholm equation of the Second 
Kind 


E‘A*,y) = + jJ c + (» - y'l) 2 ]*' (9) 

This may be solved numerically by converting the integral operator into a 
matrix operator via MoM to form the matrix equation 


U = % ( 10 ) 

where each entry of J 2 and E' t corresponds to a spatial sample. If we look at 
the impedance matrix, Z, typically it is dense. This is shown in figure 4 which 
is the impedance matrix generated by pulse basis-point matching on a lAxlA 
corner reflector where clearly all elements of this matrix are important. 

If however we apply a wavelet transform to both sides of (10), we get the 
transformed equation in the scale-domain 


zj, = c; (ii) 

Figure 5 illustrates the matrix shown in figure 4 after wavelet transformation 
by the DAUB 10 analyzing wavelet. Clearly, this matrix has fewer signifi- 
cant entries and is thus a sparse matrix. For images, edges tend to be very 
important thus we would use very compact (low-order) wavelets in order pre- 
serve these edges. However, when one has a smooth matrix such as the one 
shown in figure 4, smooth wavelets (large-order) will tend to have greater 
compaction capability than a small-order wavelet. In order to verify this 
claim, we looked at the compression realized by using different wavelets from 
the Daubechies family and a truncation tolerance of 10 -2 . Table 1 summa- 
rizes the results of this study for the comer reflector when the matrix size 
was 32x32. We observe that the apparent storage was reduced to 0(2ON) 
from N 2 without loss of fidelity as shown in figure 6. However, the large- 
order wavelets did not exhibit a significant advantage with respect to the 
small-order wavelets. 
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The reason for this unexpected result is the fact that the matrix in ques- 
tion is not large. A similar comparison is shown in table 2 for the case 
of a 6A curved(inclusion angle = 45°) metallic strip which used a I2bxl^ 
impedance matrix. We recognize some advantage using large-order wavelets 
while maintaining fidelity as seen in figure 7. As a final example we look 
at the simulation of a metallic half plane-linearly tapered resistive half plane 
junction which used a 1024x1024 matrix. A clear advantage using iarge-order 
wavelets is seen in table 3 where the apparent storage is reduced to G( - ON) 
which results in a 93% storage savings without loss of fidelity as shown m 

figure 8! ... , 

The astute engineer will recognize the use of the term ‘apparent in the 

preceding storage analysis. This is dne to the (act that yon need to generate 
the entire impedance matrix before transformation which results in 0(A ) 
storage. However, we have devised a simple block-by-block technique which 
will result in an actual memory savings. If one writes (10) as an augmented 
matrix equation where there are (JV/A') 2 blocks of A' 2 elements each where 
K is of course a power of two(a = 2) then we may write (11) as 


'£,Z,kJk = Si / = l-N/K ( 12 ) 

k-\ 

To illustrate this technique, we looked at the junction simulation previously 
given. Table 4 compares the compression results realized when A = 64 and 
the tolerance was set at 10' 4 . The total storage requirement is K 2 + 0(240 A ) 
which is a 75% storage savings compared to the original dense matrix. The 
degraded performance is due to that the fact that each block will require a 
large smooth component as well as detail components which are a result of 
the artificial boundaries created by the block method. It should be stressed 
that in spite of these artificial high components, fidelity is maintained as 
shown in figure 9. These observations are shown by comparing figure 10 
and 11 which is the impedance matrix of a 4 A flat strip when the DAUB 10 
wavelet transform is applied on the entire matrix and in a block-by- block 
fashion with sixteen 32x32 blocks, respectively. 
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4 Integral Equations 

The matrix compression technique described in the previous section is a 
rather indirect use of wavelets in electromagnetics. In this section, we will 
look at a more direct application of wavelets for the solution of integral 
equations. This method will not prove to be very useful but will provide a 
reference for future research. 

The E- polarized scattering by a metallic strip may be solved by deter- 
mining the current via the integral equation 

£z( x >J/) = J c i y )H^\k 0 \yJ(x — x') 2 + (y - j/') 2 |]d/ (13) 

where the integration domain, C, suggests the parameterization 

x = x{t) ( 14 ) 

y = y{ 0 

Thus ( 13) may be written 

£;(,) = 7 Jc. (is) 

and 

r (<>0 = V (*( 0 -*( 0) 2 + ( y ( 0 - y ( 0) 2 ( 16 ) 

Let us now employ the method of weighted residuals using the discrete 
wavelet as an entire domain weight 

X. - mT ’) dt = 7 X, X, WM - m2-‘)dt(n) 

If we now expand the unknown current in terms of wavelet basis(8) and 
realize that all of the resulting integrals are DWTs, we arrive at the matrix 
equation 


= % (18) 

where now the impedance matrix is given directly by a two-dimensional DWT 
of the Hankel function of (17) 

2 = MW[*.|r(M')||] (19) 
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and the unknown vector is composed of modal coefficients. This solution 
is analogous to finding a Fourier series expansion solution if the complex 
exponential is used rather than the wavelet. We also note that the scattered 
field is computed efficiently as 


A'-i 




7 rk 0 p 




n=0 


-j-S' 2 (n)J(n) 

4 


( 20 ) 


where the excitation is computed via the FW'T 

£•< _ yy^| c -»*o[eo»*^(0 + »* n ^v(0]j (21) 

Figure 12 illustrates that the aforementioned formulation does recover the 
solution for a 1A metallic strip when 64 coefficients are determined. This 
is an unfortunately large number of terms required. The usefulness of this 
method does not improve with problem size as seen in figure 13 where even 
1024 coefficients are not sufficient to characterize a 150A strip. The failings 
of this particular formulation will be discussed in the next section. 


5 Conclusions and Future Work 

This paper has presented several different uses for wavelets or families of 
basis functions for computational electromagnetics. The relevant features 
and characteristics of wavelets were briefly discussed. Since the scope of 
this paper precludes a detailed review of wavelets, we mention that Strang[9] 
and Daubechies[10] give excellent reviews of wavelets. Mallatfl] not only 
introduces a particularly useful interpretation of wavelets and the wavelet 
transform, he developed from this viewpoint a quite efficient Fast Wavelet 
Transform which rivals the FFT in speed. Mallat also served to present this 
field in terms which are more accessible to engineers as compared the 
mathematicians “flavor” given by Daubechies. 

We developed a novel technique adapted from image processing for gener- 
ating a sparse matrix where previously a dense matrix was required. It should 
be stressed that the utility of this method solely is based of the scale charac- 
teristics of the impedance matrix rather than some restrictive mathematical 
operation such as convolution. The block-by-block compression method pro- 
posed resulted in significant storage savings. In the future, we are interested 
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in improving this method with respect to the degraded performance caused 
by block processing. It is reasonable to anticipate some improvement since 
only the most readily implemented techniques and analysing wavelets were 
employed for this study and we understand the cause of the degradation. 
However it should be pointed out that if we can approach the efficiency of 
the entire matrix method, the block by block technique could be used with 
any existing code and a sparse matrix solver to dramatically extend the size of 
problems which could be handled with integral equation methods regardless 
of the geometry shape. 

The second method looked at was a direct solution of the integral equa- 
tion by wavelet expansion. This proved very disappointing. It is known 
that a Fourier series expansion of the current requires the retention of many 
terms. However, it was felt that the compact support property of wavelets 
would result in a dramatic reduction in the number of required coefficients. 
Although this initial trial did not succeed, we note that the wavelets cho- 
sen were available rather than the best possible. Further study may reveal 
more appropriate analyzing wavelets for our purpose or perhaps higher order 
wavelets are required for sufficient approximation. 

Another avenue of research is employing w-avelet testing with sub-domain 
basis functions. One may think of this operation as row-only compression 
rather than row and column compression which would result from the two- 
dimensional processing as described above. However, there would be no 
artificial components along the rows and by inspecting figure 11 we note 
that a significant amount of the spurious components are along rows of the 
matrix. This should result in a sparse matrix due to the local support of the 
basis and the compact support of the weights and a particularly appealing 
choice of basis are the ones used by Canningfll] since this basis will aid in 
compression along columns of the resulting matrix. To this date, we have not 
investigated this method. In addition, further study of Alpert et a/[7] may 
yield further insight. It is noted that [7] investigated non-oscillating kernel 
functions while our kernel is certainly still a matter of research. 

Wavelet processing and direct solution offers a promising avenue for ex- 
tending integral equation techniques to much larger structures. It is noted 
that all methods described herein used FWT to evaluate all the integrals in 
an efficient manner. Although the examples presented involved the MoM, we 
are interested in using these techniques when terminating a Finite Element 
mesh with a boundary integral. An efficient solution in terms of unknowns 
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may be found for arbitrarily closures. This is in contrast the few available 
closures if one employs the CGFFT method to solve the boundary integral. 
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Table 1 


lambda Cornet 

Reflector 

using 32 Unknowns 

and tol=10 / '-2 


Wavelel 

Zeros 

Non-zeros 

Total 

Order 

DAUB4 

370 

654 

1 024 

20N 

DAUBG 

451 

573 

1 024 

20N 

DAUB8 

454 

570 

1024 

20N 

DAUB10 

444 

580 

1024 

20N 

DAUB12 

439 

585 

1024 

20N 

DAUB14 

436 

588 

1024 

20N 

DAUB16 

426 

598 

1024 

20N 

DAUB18 

420 

604 

1 024 

20N 

DAUB20 

387 

637 

1 024 

20N 


Table 2 


lambda Curved 

Strip using 

128 unknowns 

and tol=10 A -2 


Wavelel 

Zeros 

Non-zeros 

Total 

Order 

DAUB4 

10311 

6073 

16384 

50N 

DAUB6 

12312 

4072 

16384 

30 N 

DAUB8 

13214 

3170 

1 6384 

20N 

DAUB10 

13388 

2996 

16384 

20N 

DAUB12 

1 3597 

2787 

16384 

20N 

DAUB14 

13591 

2793 

16384 

20N 

DAUB16 

1 3643 

2741 

1 6384 

20N 

DAUB18 

13713 

2671 

16384 

20N 

DAUB20 

1 3656 

2728 

1 6384 

20N 


Resistive Junction using 
Wavelet Zeros 

Table 3 

1024 unknowns 
Non-zeros 

and tol=10 A -4 
T olal 

Order 

DAUB4 

573834 

474742 

1048576 

460N 

DAUB6 

777806 

270770 

1048576 

260N 

DAUB8 

871382 

177194 

1048576 

1 70N 

DAUB 10 

917716 

130860 

1.048576 

1 30N 

DAUB 12 

940642 

107934 

1048576 

1 10N 

DAUB14 

954802 

93774 

1048576 

90N 

DAUB16 

966321 

82255 

1048576 

80N 

DAUB18 

973085 

75491 

1048576 

70N 

DAUB20 

978972 

69604 

1048576 

70N 


Resistive Junction using 
Wavelet Zeros 

Tabel 4 

block method w/ 
Non-zeros 

block size=64x64 
Total 

Order 

DAUB4 

268401 

780175 

1048576 

760N 

DAUB6 

432037 

616539 

1048576 

600N 

DAUB8 

569457 

479119 

1048576 

470N 

DAUB10 

685377 

363199 

1048576 

350N 

DAUB 12 

744733 

303843 

1048576 

300N 

DAUB14 

786832 

261744 

1048576 

260N 

DAUB16 

796331 

252245 

1048576 

250N 

DAUB10 

797724 

250852 

1048576 

240N 

DAUB20 

799690 

248886 

1048576 

240N 
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Figure 3. a) h 0 ( x) — X| - i/a. i/ 2 }( *)• 

b) /ij(x) - 2lw(n)X|. / 2-i/4 -/7*i/-‘l (jr) ' 

c) is decomposed into its “components", each component is a multiple of the characteristic 
function of an interval of length (The “components" of h j would have width 2~') 

d) Each “component” is replaced by a proportional version of > 1 ,. centered around the same 
point a j, the component, and scaled down by a factor {This scaling factor would be 2~ ' for h, , 

e) The functions in d) are added to constitute h 2 (or if one starts from h f in c>) 



1/vXl a. L-shaped Comer Reflector Impedance Matrix before Wavelet Transform 
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lhlA. L- shaped Comer Reflector Impedance Matrix after DAUB 10 Wavelet Transform 
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4a Strip Impedance Matrix after Wavelet Transform performed all at once 
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